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Abstract 

The capacity to predict and control bioprocesses is perhaps one of the most important objectives of 
biotechnology. Computational simulation is an established methodology for the design and optimization 
of bioprocesses, where the finite elements method (FEM) is at the state-of-art engineering multi-physics 
simulation system, with tools such as Finite Element Analysis (FEA) and Computational Fluid Dynamics 
(CFD). 

Although FEA and CFD are currently applied to bioreactor design, most simulations are restricted 
to the multi-physics capabilities of the existing sofware packages. This manuscript is a contribution 
for the consolidation of FEM in computational biotechnology, by presenting a comprehensive review of 
finite element procedures of the most common enzymatic mechanisms found in biotechnological processes, 
such as, enzyme activation, Michaelis Menten, competitive inhibition, non-competitive inhibition, anti¬ 
competitive inhibition, competition by substrate, sequential random mechanism, ping-pong bi-bi and 
Theorel-Chance. 

Most importantly, the manuscript opens the possibility for the use of FEM in conjunction with in-silico 
models of metabolic networks, as well as, chemical networks in order to simulate complex bioprocesses in 
biotechnology, putting emphasis into flux balance analysis, pheno-metabolomics space exploration in time 
and space, overcoming the limitations of assuming chemostat conditions in systems biology computations. 
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Introduction 

Predicting the behavior of bioprocesses is one of the major goals of biotechnology. Computational simulation 
is today a valuated tool for predicting, monitoring and controlling the status of fermentations, as well as, for 
optimizing fermentation conditions, minimizing trial and error experimental procedures. 

Computational design is recognized as a standard prototyping tool outside the bioengineering area (e.g. 
automotive and aviation), where it significant reduces costs during design, prototyping and testing phases. All of 
these, generally involve high experimental load and trained personnel in different areas of research and engineering. 
The same is also becomming a reality in biotechnology with the advent of systems and synthetic biology. 

Traditional experimental methods are limited by the number of recorded parameters for a holistic systems 
characterization. The conjunction of high-throughput methods (e.g. mass spectroscopy, microarrays, sequencing, 
spectroscopy and electrochemistry) are today elected for validation of state-of-the-art ’in-silico’ chemical and 
genome scale models (GSM). Computational simulation provides detailed information in time and space. The 
Finite Element Method (FEM) is of the ’’heart” of many Finite Element Analysis (FEA) and Computational Fluid 
Dynamics (CFD) software for simulating physical phenomena, such as, heat transfer, mass transfer, radiation, 
fluid dynamics, structural and elasticity, but it can also be used in biotechnology for simulation of chemical, 
biochemical reactions, and cellular dynamics m- 

FEM was not initially developed for computational biology and bioprocesses simulation. It has been devoted 
to industrial prototyping of biotech machinery [3H5]. It is not yet usual the application of FEM for the simulation 
of complex biological or chemical systems [SHS]. In sophisticated developments, FEM has been used to compute 
microscopic properties, such as: i) the study of membrane elasticity [10]: ii) electrostatic interactions between 
proteins m-, iii) mechanical modeling of ion channels m-, iv) applying FEM in microscopy for physical properties 
estimation m- FEM has also been applied to the study of enzyme kinetics by continuous diffusional biomolecular 
systems given by the Smoluchwski equation. It has proven to be a good alternative to the traditional spherical 
criterion model, by allowing to study the complex enzyme geometries |14m8 |. 

The continuity of FE facilitates the inclusion of other phenomena such as, fluid flow, heat/mass transfer, 
electromagnetic field, forces and elasticity. The computational cost is less when computing large scale problems 
described by differential equations, where continuous solutions are common in physical phenomena even at small 
scales (e.g. force fields, diffusion, heat transfer) | 19H22| . 

The main steps in FEA involve: i) Pre-processing; ii) Resolving the PDEs or ODEs in the physical-time 
domain; and iii) post-processing. Pre-processing generally involves: i) ensure that PDEs and ODEs are interactive 
for multi-physics and chemical, biochemical and microbiological models; ii) ensure that the solution is stable and 
accurate in the physical-time domain by optimizing the mesh refinement and time steps from computer assisted 
design software |23ff25| or in more complex geometries (e.g. biological tissues) by digital scanning and 2D/3D 
reconstruction methods. This methodology has a number of advantages, such as the treatment of problems 
on complex irregular shapes, non-uniform meshing to reflect different levels of multi-scale detail, treatment of 
boundary conditions using continuous solutions and the construction of higher-order approximations to improve 
accuracy of numerical solutions. Both biological materials, as well as, bioreators display irregular geometries and 
non-homogeneous physical-chemical properties, which makes difficult to sustain a chemostat hypothesis. FEM 
not only overcomes such hurdle, but when used in conjunction with inverse problems makes possible to minimizing 
the error between simulation and experimental datasets obtained in discrete positions of space, to improve model 
predictions |26H28| . As biological processes implie multi-physics and multi-scale simulations, it becomes essential 
to: i) develop the correct relationship between physical-chemical, biochemical and microbiological models; ii) 
ensure that all used model parameters are correctly determined against experimental data by inverse methods 
and statistical analysis |29| . 

High-throughput molecular biology and analytical chemistry technologies are exponentially increasing chemical 
and biological ’omics’ information databases (e.g. genomics, metabolomics, transcriptomics, proteomics and 
protein interactions) (see Figure [T]). The available information allowed the emergence of the annotation of gene, 
protein and metabolic functions, as well as, regulatory mechanisms, so that, network reconstructions of complex 
biological systems are today feasible. Network models gave rise to the development of ’in-silico’ network organisms, 
reconstructed from curing the information present in both databases and publications |30I31| . allowing the analysis 
of network properties and topology, as well as, the comprehensive analysis of cellular functions by systems biology 
approaches |32| . 

Connecting all mathematical models on a multi-scale and multi-physics strategy is one of the most important 
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challenges for understanding the complexity of chemical and biological systems [^. This manuscript is a con¬ 
tribution for the basis of the use of the finite element method procedures for the integration of enzyme kinetics, 
chemical and genome scale network models (’in-silico’ strains) as a complex systems multi-scale and multi-physics 
computational modeling research area. This communication is not a comprehensive presentation of the finite 
elements method, and therefore background on numerical modeling is necessary to make use of the presented 
equations. 


Materials and Methods 

The finite elements method 

The FEM is considerably different from the most common discretization methodologies, such as Finite Differences 
(FD), Finite Volumes (FV) and Lattice-Boltzmann (LB) methods. Although elements are geometrically equal, 
FEM ensures that the solution is continuous inside each element, solved by a weak solutions to a variational 
optimization of a quadratic problem, being the solution inside a physical given by a piecewise continuity - the 
shape function. 

The following steps resume the FEM methodology: i) Passing from global to local coordinates for the shape 
function of such element; ii) Variational analysis - determining the solution to the variational problem by weakening 
the solution inside the finite element; iii) matrix assembly of all equations; and iv) solving |26il27l[55U38| and 
rendering results into graphical mode [281I371I39] . 

The variational method 

Changing a Partial Differential Equation (PDE) or an Ordinary Differential Equation (ODE) into the variational 
form is the main procedure for any FEM discretization. The simplest form of a variational (V(a:)), states for two 
continuous functions h{x) and v(x)'. 


V{x)= f h{x)v{x)dx = 0 (1) 

J a 

which means that v(x), a weighting or testing function, can be chosen to force the residuals h{x) to be 
zero inside the finite element interval [a,ti]. The variational problem is posed in the finite element space (D). 
The variational can be solved by the direct substitution of the residuals function {h{x)) and weighting function 
(u(a;)) and minimization (Garlekin’s method) or by the minimization of a linear functional (functional method) 
|34[|4UfH2| . The variational method states that there is a solution to the problem of eq. [T] given by: 

B{u,v) = £{v) (2) 


which satisfies the solution u = u*, for any trial function v (or shape function). B{u, v) is a bi-linear functional 
dependent upon the original and trial function and C{v) a linear functional dependent only of the trial function. 
The condition above is only possible to be obtained if the following functional is minimized in the case of 1st 
order differential equations: 


V(u, u) 


^B(u,u) - £{u) 


(3) 


In order to u (^ = 0) using the shape function {u = Nu), holds the solution to B{u, v) — £{v), and to the 
variational problem in eq. [T] [34l|40l|4T] . 

Shape Functions and Elements 

Shape function is the continuous approximate solution to the variational problem inside the finite element space. 
The shape function is only dependent upon the type and shape of the finite element. Elements can be grouped 
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into their different interpolation functions: i) first order - linear elements; ii) second order - quadratic elements, 
and iii) third order - cubic elements (higher order shape functions are unusual). 

The most common FE shape functions are presented in Table [T] These describe a piecewise solution to the 
variational problem in the physical domain. The FEM Is generally used in the natural coordinate system [33]. 
For example, the subtract concentration Cs inside a linear rectangular element: 


Cs{x, y) = NiCs^i + NjCs,j + NkCs,k + NmCs,m (4) 

where Cs{x,y) is the enzyme concentration in each of the physical positions inside the element, and Cs,i, 
Cs.j, Cs,k, Cs,m at the nodal presents i, j, k and m respectively, Ni, Nj, Nk and Nm are the shape coefhclents 
and are dependent on the elements coordinates ElllSlIlI]. 

The shape function can be presented in the compact matrix format: 


CE{x,y)=[Ni Nj Nk Nrr.]- 


CE,i 
Cej 
CE,k 
CE ,m 


= N-Ce 


( 5 ) 


where Ni ... Nm are the shape coefficients and Ce.i ■ ■ ■ CE,m the enzyme concentrations at the element nodes 
i to m, respectively. The discretization presented in this manuscript can be further extended to the different types 
of elements using similar mathematical reasoning. 


Reaction Models 

Reactions in space-time 

Diffusion dependent enzyme reactions are well described by the 2nd Fick law: 

^-V(DVC7)-r* = / (6) 

where C is the specimen concentration {mole.dm~^), D the mass diffusivity f the force vector and 

r* the reaction rate of C per unit value {mole.s ~^The simplest reaction term r* to be added to eq[6]is 
the first order kinetic: 

A^B 


which can be described by: 


^ = -kCa (7) 

where the kinetic rate k (s“^) is a function of temperature given the modified Arrhenius law. Such states that 
the decay of A is proportional to the probability of finding A {p{A)) molecule inside the finite element space, that 
is p{A) oc Ca- Consequently, the first order reaction is proportional to Its concentration inside the finite element. 
For the sake of simplicity, lets assume this reaction is occurring inside a linear triangle where C'(A) Is given by: 

Ca = NiCa,i + ^jCa,j + NkCa,k ( 8 ) 

which is the space distribution of probabilities of finding (A) inside the finite element space. Therefore, Ni, 
Nj and Nk map the random movements of A molecules inside the finite elements, proportional to speciemens 
concentration. 

Once formulated the variational problem is possible to obtain: 

Vn = J + kCa^ ■ v{n)dn 

= [ ^ ■dD.+ [ kCa- v{n)dD. = 0 
Jn Jn 


(9) 
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Under these circumstances the variational can be solved by using a linear functional which holds the true 
solution after the minimization of the bilinear functional V{Ca,v)'. 

V{Ca,v)^^B{Ca,v)+C{v) (10) 

where the functionals B{Ca,v) and £(«) are given by: 

B{Ca,v)= [ kCa-v{Q.)d9. 

Jn 

c{v) = / ^ • v{Q.)dn 
Jn 

and therefore the functional V[Ca) takes the form of: 

v{Ca,v) = J ^kCadn + 1 (^^yCadn 

Substituting the element functions in the first term of eq. 1121 yields: 


V — o / ^dSliki + Njkj + Nkkk') • 

^ Jn 

{NiCa,i + + dJkCa,k) dll (13) 

That once minimised for the node i, holds: 


( 11 ) 

( 12 ) 


= [ [{hN! + kjNfNj + kkNfN^) ■ Ca,i 
OUa,i Jn 

+ {kiNfNj + kjNiNf + kuNiNjNk) ■ Ca,j 
+ {hNfNk + kjNiNjNk + kkNiNl) ■ CaM (14) 

The same minimization is necessary to be made in terms of Ca,j and Ca,k, to obtain all the elements of the 
final matrix K. After algebraic manipulation, the stiffness matrix (K) is possible to be described in the matrix 
format by: 

K = f (15) 

Jn 


where k is the column vector k = [ki kj kk], and the kinetic rate inside the finite element is given by Nk. If 
one considers fc as a row vector, than the solution is K = N^Nk*NcAl, since these are symmetric matrices. 
Similarly, for the second term of eq. 1121 


V = 



dCa,i 

dt 


+ N, 


dCa,j 

dt 


+ Nk 


dCa,k \ 

dt J 


{NiCa,i + NjCaJ + NkCa,k) 80 


That once minimised in terms of Ca,i yields: 


av 

dCa.i 



+ NiNj 


dCa,j 

dt 


+ NiNk 


dCa,k \ 

dt J 


80 


(16) 


(17) 
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The same kind of minimization is necessary for Ca,j and Ca,k, to obtain the final matrix that will enable the 
FEM method computation. After algebraic manipulation, the full minimization of the variational is possible to 
be presented in the matrix format: 


dV f 

— = / • Ca (18) 

oCa Jq 

Therefore, the chemical reaction can be computed accross the physical domain by: 

[ • Ca + / Kan • Ca = 0 (19) 

Jn Jq 

where N'^NdQ presents the probabilities of random movements of the molecule a in any direction inside the 
finite element and Kdfi the probabilities of effective conformational changes of a across the finite element. An 
important assumption in this discretization, is the fact that kinetic rate is not constant across the physical domain. 
Such occurs in non-homogeneous biological materials. If one considers constant kinetics, than K = kc N^NSfl, 
where ko is a constant kinetic rate. For more chemical reaction mechanisms, please consult [5]. 

Although the 1st order reaction kinetics is the most simplest mechanism, it is still the most widely used to 
represent both systems chemistry and ’in-silico’ organisms, where single steps are considered uni-molecular, and 
in the last case, catalyzed by an enzyme, being possible to be used in conjunction with reaction networks and 
GSM. 


Second order kinetics 

The simplest form of reaction given by molecular colisions, is the second order reaction kinetic: 

A + B^C 

where, 

^ = -kCaCb (20) 


Ca and Cb are the concentrations of A and B specimens inside the finite element. In this case, reaction 
only occurs once there are effective collisions between A and B. Therefore, inside any linear finite element the 
variational form is presented as follows: 


V = [ '^v{Q.)d^ f kCaCbv{fl)dQ 

Jq dt 2 Jq 

where k, Ca and Cb vary consistently inside the finite element, taking the form: 


( 21 ) 


ViCa) = [ ^CadQ -b i / kCbCldn 

Jq dt A Jq 

which minimizing for node i, j and fc, attains: 

^ r f kcCaNidft 

oCa,i Jn dt Jn 

Which after the variational minimization, the solution yields: 


( 22 ) 


(23) 


[ N'^Ndn ■Ca+ [ Mdn ■ CaCl ■ u = 0 (24) 

Jn Jn 

where, M = diag(N)K, and U is the column vector [1 1 1]. M expresses the frequency of A and B to 
react inside the finite element. Afterwards, both equations for Ca and Cb solution must be computed with both 
equations. Moreover, the term CaCb expresses all possible collision probabilities between a and b specimens 
inside the finite element. The same is possible to derive for auto-catalyzed reactions {A-\-A^C), being possible 
to show that the solution is held by: N^Ndfl ■ Ca + Mdn • CaC^ ■ U = 0 . 
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Results and Discussion 


Enzymatic models 

Enzyme activation 

The enzymatic activation/inactivation is an example of fractional conversion model [43], that describes an equi¬ 
librium between two species and E, which correspond to inactive and active enzymes, respectively. 

E° ^E 

In this case, the concentration of E° and E is established by a dynamical equilibrium by: 


-h feiCgo - k-iCE = 0 
at 

^ -h = 0 


(25) 

(26) 


After manipulations, the following finite element formulation inside the linear finite element for inactive and 
active enzymes, respectively: 


• Ceo + / -Ceo- / Ei-idn • Cb = 0 


Reaction chain 

Consider the following reaction chain: 


Ai. 


Ai —^ Ai 


+1 


where i is the i’th specimen in the reaction chain. By direct comparison with previous formulations it is simple 
to derive the FEM formulation for each specimen: 


N^Nan-Ca;- / KiSfl-Cai-f 


/ 


idQ ■ Cui-i = 0 


The full reaction chain is computed by joining all the specimens equation matrixes. 


Michaelis-Menten model 

The most widely known enzymatic model is the Michaelis-Menten mechanism: 


E + S^ES^E + P 

which can be expressed by a balance to each species: 

ACes 


dt 


— UxCeCs -I- k-iCES + k2CES = 0 


dCs 

dt 


-I- kiCsCs — k-iCES = 0 


dCp 

dt 


— k2CES = 0 


(27) 

(28) 
(29) 


ACe 

dt 


-I- klCECs — k-lCES — k2CES = 0 


which inside the finite element can be expressed as: 


(30) 






• Ces - / Mian • Ces ■ Ci ■ U 


+ / K-iSn • Ces + / KaSn • Ces = 0 

Jq Jn 


/ N-' Ndn • Cs + / Mdn • Ce ■ ■ U 

Ifi Jci 

- f K-idn-c^s = 0 

Jn 


[ N^Ndfl-Cp- [ K2a^l-CES = 0 
Jn Jn 
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(31) 


(32) 


N^Nan • Ce + / Man • Ce ■ cj • u 


- f K_i Ces - [ Kaan • Ces = 0 

Jn Jn 


(33) 


Inhibition of enzymatic activity 

Enzyme inactivation is both a control mechanisms, as well as, a lowering yield factor by an inhibitor I (natural 
or synthetic) which interacts with the enzyme, decreasing the catalytic activity. 


Competitive inhibition 

Competitive inhibition occurs when an inhibitor (/) (Figure [2| competes with the substrate for the active center, 
being represented by: 

E + S^ES^E + P 
E + I ^ El 


which can be expressed by a balance to each species: 

ACes 


dt 


— kiCECs + k-iCss + k2CES = 0 


dCs 

dt 


+ kiCECs — k-iCss = 0 


dCp 

dt 


— k2CES = 0 


dCs 

dt 


+ kiCsCs — k-iCES — k2CES + k^CECi — k-zCEi = 0 


dCsi 

dt 

dl 


— kiCECi + k-iCEi = 0 


dt 


+ kiCECi — k-iCEi = 0 


(34) 

(35) 

(36) 

(37) 

(38) 

(39) 


and after solving the variational problem, the solution can be expressed as: 




9 


/ N^N9n-CEs 

Jn 

- I Mi9n ■ Ce ■ Cs • U 
Jn 

+ [ K^idn ■ Ces + [ K2dn ■ Ces = 0 
Jn Jn 


N^NdO. • Cs + / Mi(9fi • Ce ■ Cs • U 


- / K^ian cEs = 0 


[ Cp- [ Kaan • Ces = 0 

jQ, Jn 


(40) 


(41) 


(42) 


N^Nan ■ Ce + / MiSn ■ Ce ■ Cj ■ U 


- [ K-idn Ces - / KaSH • Ces 
Jn Jn 


+ / MaSn • Ce • Ci" • U - / K-gail • Cei = 0 
Jn Jn 


N^Nan - Cei - / MaaH • Ce • C?' ■ U 


+ / K_3an • Cei = 0 
Jn 


(43) 


(44) 


N'^Nan • Cl + / Mzan • Ce • c?' • U 


- / K_3afl ■ Cei = 0 
Jn 


(45) 


where Mi and M 2 express the frequency of E-S and E-I to react inside the finite element. 


Non-competitive inhibition 

Non-competitive inhibition occurs when an inhibitor (7) reversibly establishes a chemical bound with the enzyme 
which is not the active site, but nevertheless affects its catalytic activity, being possible to be expressed by the 
mechanism: 


E + S^ES^E + P 
E + I ^ El 
ES + I^ ESI 


EI + S^ ESI 
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which can be expressed by a balance to each species; 
dCss 


dt 


— klCsCs + k-lCsS + k2CES + kiCEsCl — k-4CESI = 0 


dCs 

dt 


+ kiCECs — k-iCES — k-^CESi + k^CEiCs = 0 


dCp 

dt 


— k2CES = 0 


dCE 

dt 


+ kiCECOs — k-iCss — k2CES + k^CECi — k-sCEi = 0 


dCEi 

dt 

dCssi 

dt 

dCi 

dt 


— kiCECi + k-zCEiCs + ksCEiCs — k-^Cssi = 0 
+ k-sCESi — k^CEiCs + k-iCESi — kiCEsCi = 0 
— k-ACESi + kACEsCi — k-sCsi + ksCsCi = 0 


which inside the finite element can be expressed as: 


■ Ces - / Mian • Ce ■ cl • U 


+ / K-iSn • Ces + / K2an • Ces 

Jn Jn 

+ [ Mzan • Ces • C| • U - [ K-aO^I ■ Cesi = 0 
Jn Jn 


/ Nan • Cs + / Mian • Ce • c| ■ u 

In Jn 

- [ K_ian ■ Ces - / K.san ■ Cesi 

Jn Jn 


Maan • Cei ■ C| • U = 0 


/ N-* Nan • Cp - / Kaan ■ Ces = o 
In Jn 


+ [ Mgan • Ce • c| ■ u - / K_3an • Cei = o 
Jq Jq 


[ N^Nan • Ce + / Mian • Ce • c| • u 
Jn Jn 


- / K_i • Ces - / Kaan • Ces 

Jn Jn 


N^Nan - Cei - / Mian • Ce ■ C| ■ U 


+ / Maan- Cei ■ c| • u + / M4an • Cei • c|u 


- / K_5an-CEsi = o 


(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 


(53) 


(54) 

(55) 


(56) 


(57) 
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/ N-* Nan • Cesi + / K_6an-CEsi 
in in 

[ M 4 an • Cei ■ cj • u + [ K-4an-CEsi 

In Jn 

- [ M4an ■ Ces ■ c?' ■ u = 0 
Jn 


(58) 


/ N" Nan • Cl - / K_ 4 an • Cesi 
in in 

+ [ M4an • Ces • C? • U - [ K-adQ. ■ Cei 
Jn Jn 

+ [ Mgan • Ce • c?’ • u = 0 

Jn 


(59) 


where Mi, M 2 , M 3 , M 4 express colliding probabilities of — 5, ES — I, E — I and El — S. 

Anti-competitive inhibition 

When an inhibitor links itself reversibly to enzyme-substrate complex and not to the free enzyme, this is known 
as anti-competitive inhibition mechanism: 

E + S^ES^E + P 
ES + I^ ESI 

which can be expressed by a balance to each species: 
dCss 


dt 


— kiCsCs + k-iCss + k2CES + kaCssCi — k-gCESi = 0 


dCs 

dt 


+ kiCECs — k-iCES = 0 


dCp 

dt 


— k2CES = 0 


dCE 

dt 


+ kiCECs — k-iCES — k2CES = 0 


ACesi 

dt 

dCi 


— kaCEsCi + k-aCESi = 0 


dt 


+ kaCEsCi — k-aCESi = 0 


which inside the finite element can be expressed as: 


N^Nan ■ Ces - / Miatl • Ce ■ C| • U 


(60) 

(61) 

(62) 

(63) 

(64) 

(65) 


+ / K-iail • Ces + / ■ Ces 

Jn Jn 

+ [ M2dn ■ Ces • C?' • U - [ K^sdn ■ Cesi = 0 
Jn Jn 


(66) 
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• Cs + / Mian • ci • Cs • u 


- / K_ian ■ Ces = 0 


(67) 


/ N" Nan • Cp - / Kaan • Ces = o 

Iq Jq, 


( 68 ) 


/ N'^Nan ■ Ce + / Maan ■ Ce ■ ci' • u 

/n 

- [ K-idn ■ Ces - / • Ces = 0 

Jq Jq 


(69) 


/ N-' Nan • Cesi - / Mian • Ces • ci • u 
/ n Jn 

+ f K_3an • Cesi = 0 


(70) 


N'^Nan ■ Ci + / Maan ■ Ces ■ Ci • U 


- / K-aan • Cesi = 0 


(71) 


where Mi and Ma express the frequency of E — S and ES — I colisions the finite element. 

Ping-Pong Bi-Bi mechanism 

In Ping-Pong Bi-Bi mechanisms, one of the substrates connects to the enzyme and one of the resulting products 
releases before the second substrate can connect: 

E + EA^ E*P ^ E* +P 

E* +B ^ E*B ^ EQ ^ E + Q 

which can be expressed by a balance to each species: 


N'^Nan • Ce + / Mian ■ Ce • cX ■ u 


- / Kean • Ceq = 0 


(72) 


N^Nan ■ Ca + / Maan • Ce • ci • u = o 


(73) 


N^Nan • Cea - / Mian • Ce ■ cX • u 


- / K_2an • Ce*p + / Kaan • Cea = o 

Jq Jq 


(74) 
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/ N-' N9fl • Ce* - / KadQ ■ Ce*p 

Jn Jn 

+ [ Msan ■ Ce* • Ci ■ U = 0 


/ N-'Nan-CE*p- / K_2an-CE*p 
In Jn 

+ f Kadn ■ Ce*p = 0 


/ Nan • Cp - / KaSn • Ce*p = 0 
In Jn 


N'^Nan ■ Ce + / Mian ■ Ce • cl ■ u 

! in 


- / Kean • Ceq = 0 


N'^Nan • Ca + / Mzan • Ce ■ Cs • u = o 


/ N^Nan • Cea - / Mian ■ Ce ■ cj ■ u 

/n in 

- [ K-2dn • Ce^p + [ K2dn • Cea = 0 

Jq Jq 


/ N^NdQ-CE* - Kad^ CE*^ 

Iq Jn 

+ [ Msdfl • Ce* • Ci • U = 0 

Jn 


/ N-'Nan-CE*p- / K_2an-CE*p 
/n in 


+ / Kaan • Ce*p = 0 


[ N^Nan • Cp - / Kgan • Ce*p = o 

Jn 


N'^Nan • Cb + / Msan • Ce* • cl ■ u = o 


(75) 

(76) 

(77) 

(78) 

(79) 

(80) 

(81) 

(82) 

(83) 

(84) 
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• Ce*b - / Mgan • Ce* • ci ■ u 


+ / K5an CE-B+ / K_59fl CEQ 

Jn Jn 

[ Kedn ■ Ceq = 0 


(85) 


/ N^'Nan • Ceq - [ Ksdil ■ Ce.b 

Jq Jq 

+ [ K-6dn • Ceq + [ KeOn • Ceq = 0 

Jn Jn 


( 86 ) 


/ N-" NSn • Cq - / • Ceq = 0 

Jn Jn 

where Mi, M 2 , M 3 expresses the collision probabilities o^ E — A, E — S, and E — respectively. 


(87) 


[ N^Nafl-CB+ / 

Jn Jn 


N'^Ndn -66+/ Mgan • Ce* ■ ci ■ u = 0 


( 88 ) 


• Ce*b - / Mgan • Ce* • ci ■ u 


+ / K5an CE*B+ / K_5afl CEQ 

Jn Jn 

[ Kedn ■ Ceq = 0 


(89) 


/ N-* Ndft ■ Ceq - / Ksdft ■ Ce.b 

Jn Jn 

+ [ K-sdil ■ Ceq + [ Kedil ■ Ceq = 0 

Jn Jn 


(90) 


[ N^Nan • Cq - / Kean • Ceq = 0 (91) 

Jn Jn 

where Mi, M 2 , M 3 expresses the collision probabilities of E — A, E — S, and E — respectively. 


N'^Naii • Ce + / Miaii ■ Ce • cX ■ u 


- / Kean • Ceq = 0 


(92) 


N'^Nan ■ Ca + / Maan • Ce • cl • U = 0 


(93) 
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• Cea - / Mian ■ Ce ■ cl ■ u 


- [ K-2an cE*p+ / K2an cEA = o 
Jn Jn 


/ N'Nan cE^-/ Kgan cE^p 

In Jn 

+ [ Mgan • Ce* • ci ■ u = 0 

Jn 


/ N"Nan-CE*p- / K_2an-CE*p 
In Jn 

+ f Ksdn ■ Ce*p = 0 


/ N" Nan • Cp - / Kgan • Ce*p = o 

In Jn 


/ Nan • Cb + / Msan • Ce* • Cb ■ u = o 

In Jn 


N^Nan • Ce*b ” / Msan ■ Ce* • ci ■ u 


+ / K5an cE*B+ / K_5an cEQ 


• Ceq = 0 


/ N' • Ceq - / • Ce*b 

In Jn 


+ / • Ceq + / • Ceq = 0 

in Jn 


(94) 


(95) 


(96) 


(97) 


(98) 


(99) 


( 100 ) 


[ N^Nan • Cq - / Kean • Ceq = o (loi) 

Jn Jn 

where Mi, M 2 , M 3 expresses the collision probabilities of -B — j 4, — S', and E — B, respectively. 

Ping-Pong Bi-Bi with parallel pathway 

In some cases, parallel pathways as in Ping-Pong Bi-Bi, such as for DD-carboxypeptidases [44], being an important 
reaction pattern to be discretized into FEM. 

E + P ^ ED ^ E*P E* +P 

E* +A^ E*A ^ ET ^ E + T 
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E* +B E*B ^ EC E + C 

the finite element formulation is presented as: 


/ • Ce + / Mian • Ce ■ eg ■ U 

In Jn 

- [ Kean • Get - / Kgan • Cec = 0 
Jn Jn 


N^Nan • Cd + / Mian ■ Ce • eg ■ u = o 


N^Nan ■ Ced — / Mian ■ Ge • eg ■ u 


- / K_2an • Ge*p + / Kaan ■ Ced = 0 

Jq Jq 


/ N"Na^^ •CE*p- / K2dn- Ced 

Jn Jn 

+ [ K-2dn CE*p+ [ Kadn-CE*^ =0 

Jn Jn 


■/: 


I N^'Nan • Ce* - / Kgan • Ce*p 
Jn Jn 

+ I Maan • Ce* • cl ■ U + [ Madil ■ Ce* • Cg • U = 0 

Jn 


/ N" Nan ■ Cp - / Kaan • Ce*p = o 

In Jn 


N^Nan • Ca + / Maan • Ce* • cl ■ u = o 


N^Nan • Ce*a - / Maan • Ce* • cj 


- / K_6an • Cet + / Ksan-CE-A = 0 
Jn Jn 


/ N-'Nan- Cet — / Ksan - Ce*a 
Jn Jn 

+ [ K^san- Cet + [ Kean- Cet = 0 

Jn Jn 


/ N" Ndn • Ct - / Kedn • Cet = 0 

In Jn 


( 102 ) 

(103) 

(104) 

(105) 

(106) 

(107) 

(108) 

(109) 

( 110 ) 

(111) 
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N^N9fl ■ Cb + / Madn • Ce* - Ci • U = 0 


( 112 ) 


• Ce*b - / MsSn ■ Ce* • ci ■ U 


- / K-sdn ■ Cec + / KadQ ■ Ceb = 0 

Jq Jq 


(113) 


/ N" Nan • Cec - / Kgan • Ceb 

Iq Jq, 


+ / • Cec - / • Cec = 0 

Jq Jq 


(114) 


/ N^Nan-Cc - / Kgan ■ Cec = o (115) 

Jn Jfi 

where Mi, M 2 , M 3 express the collisions probabilities oi E — D, E — A and E — B. 

Ternary-complex mechanisms 

Ternary-complex mechanism is also common in cellular processes (e.g. DNA polymerase). In this type of enzyme, 
two substrates need to link to the enzyme to form a ternary complex, either in sequence or random, with the 
following set of reactions: 

E + A^EA 
E + B ^ EB 
EA + B ^ EAB 
EB + A^ EAB 
EAB ^ EPQ 
EPQ EP + Q 

EPQ EQ + P 

EP ^ E + P 
EQ^E + Q 

which inside the finite element can be expressed as: 


-,T 

^B 


- / K4an ■ Ceq = 0 


u 


Iep 

(116) 

= 0 

(117) 
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Ndn • Ca + / MiOn • Ce ■ Ci - U 
! Jn 

+ j K4(9r2 • Ceb = 0 
Jn 

N'^NdD. • Cb + / Maan • Ce • Ci ■ U 

i Jn 

+ f Ksd^ • Cea = 0 


/ 


N^Ndn ■ Cea - / MiSn • Ce • Cj • U 


+ / MsSfi • Cea ■ C^ • U = 0 


/ 


Ceb — / MaSn • Ce • ci • U 


+ / M4an • Ceb ■ CI ■ U = 0 


/ 


N^Nan ■ Ceab - / Maan ■ Cea • Ci • U 


- / M4ar2 - Ceb ■ C^ ■ U+ / Ksaj^-CEAB 


/ 


- / K^san • Cepq = 0 


N' Nan • Cepq - / Ksdfl ■ Ceab 


+ / K-san • Cepq + / Kean • Cepq 


+ / Kran • Cepq = 0 


N" Nan • Cep - / Kean • Cepq 


+ / Kgan • Cep = 0 


N' Nan • Ceq - / Kran • Cepq 


+ / Kgan • Ceq = 0 


N" Nan • Cp - / Kgan • Cep = o 


/ N' Nan • Cq - / Kean • Ceq = o 

Jn Jn 

where Mi, M2, M3, M4 express the colision probabilities ot E — A, E — B, EA — B and EB- 


(118) 

(119) 

( 120 ) 

( 121 ) 

( 122 ) 

(123) 

(124) 

(125) 

(126) 
(127) 

A, respectively. 
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Rapid-equilibrium random mechanism 

In this mechanism, the enzyme is capable to randomly link to four different substrate (A, B, D) to form the 
complex EDA or EDB, producing the different molecules P, T and C [45], as follows: 

E + EA + D ^ EDA ^E + P + T 

E + D^ED + A^ EDA 

ED + B ^ EDB ^E + P + C 

E + B^EB + D^ EDB 

which inside the finite element can be expressed as: 


/ Nan ■ Ca - / K_ian ■ Cea 
In Jn 


+ / Misn ■ Ce ■ ci • u - / K4an • Ceda 

Jn Jn 


+ / M 4 an ■ Ce ■ ci ■ u = 0 

Jn 


/ N-* Nan • Ce - / Kian • Cea 

Jn Jn 

+ [ Mian • Ce ■ cl ■ U - [ KsdQ ■ Ceda 
Jn Jn 

+ [ O-sdn Ce cI -Ct- [ K4an • Ced 
Jn Jn 

+ [ M4an • Ce • cS • u - [ Kran • Cede 

Jn Jn 

+ [ O-ran ■ Ce ■ cj ■ Ct - [ Ksdn • Ceb 

Jn Jn 

+ [ Mean ■ Ce ■ ci ■ U = 0 
Jn 


/ N" NdQ • Cb - / Ksdfl • Ceb 
Jn Jn 

+ [ Msdn • Ce ■ cl • U - [ Kaan • Cede 
Jn Jn 

+ [ Mean • Ced • Ci • U = 0 
Jn 


/ N" Ndn • Cd - / • Ced 

Jn Jn 

+ [ M 4 an ■ Ce ■ cl ■ U - [ Kodn ■ Cede 
Jn Jn 

+ [ Mgan • Ceb ■ cS ■ U = 0 
Jn 


(128) 


(129) 


(130) 


(131) 
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N-' N9n • Cea - / Mian • Ce • Ci • U 


+ / K^iSn • Cea - / K-2dQ ■ Ceda 


+ / MaSn CEA-CD •U = 0 


N^Nan • Ced - / M4an • Ce • eg • u 


/ 


K_4an cED- / K-san cEDA 


+ / Mean- Ced ■ ci ■ u - / K-ean • Cede 


+ / MeOn • Ced • C^ • U = 0 


N • Ceb / K_8<9f^ ■ Ced 


/ 


+ / Msdn • Ce ■ cl • U 


/ 

Jn 


- / K^gan • Cede 


+ / Mgan • Ceb ■ Cg ■ U = 0 


N' Nan -Ceda - / Maan • Cea ■ Cg ■ U 


+ / K_ 2 an • Ceda + / Ogan • Ce ■ Ci ■ Ct 

in in 

+ f Kaan • Ceda - f Mean • Ced • CJ ■ U 


+ / K-san • Ceda = 0 


N-' Nan • Cede + / Mean • Ceb • Clg • U 


K^ean • Cede — / O -7 an • Ce • c? ■ Cc 


Kran • Cede — / Mean- Ceb • Cd • u 


+ / K—gan • Cedb = 0 


N' Nan • Cp - / Kaan ■ Ceda 


+ / O-aan • Ce • c? • Ct - / Kran • Cedb 


- / o_ran ■ Ce ■ c? ■ Cc = 0 


(132) 


(133) 


(134) 


(135) 


(136) 


(137) 
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N-'Nari-Cx- / Ka^n-CEDA 


+ [ ■ Ce ■ C? ■ Ct = 0 

Jn 


(138) 


[ N^Nan ■ Cc - / KrdQ ■ Cedb 

Jn Jq 

+ [ O^rdn • Ce • Cj ■ Cc = 0 (139) 

Jq 

where M expresses the frequency of reaction of metabolisms inside the finite element. 

Chemical networks 

When reactions are put together to describe a chemical system, it can be formalized as a graph, where reactions 
are links or edges and specimens are nodes (Figure m. Take for example the following chemical set of reactions: 

H + HCl -^H2 + CI 
HCl + 0 ^ Cl + OH 
HCl + OH ^Cl + H 2 O 

that can be represented by the graph in Figure |4] (a). In this network, all reactions involve a second order 
reaction kinetics mechanism, following the FEM discretization presented in section „ If no spacial variation is 
considered, the differential equation for the presented reaction network is as follows: 


C + S-C = 0 


(140) 


where C is the reaction rate and S the stoichiometry matrix derived from both stoichiometry and reaction 


graph. For this reaction network, the following system of equations is obtained: 




• 1 

0 

0 • 

^HCL 


1 

1 

1 



-1 

0 

0 

^Cl 

+ 

-1 

-1 

-1 

Co 


0 

1 

0 

Cqh 


0 

-1 

-1 

. CH 20 - 


0 

0 

— 1 . 


kiCn Chcl 
k^Co OhcL 
kiCoH^HCL 


( 141 ) 


which must be solved by optimization methods. The reaction network can also be represented by an incidence 
matrix (speciemens relationships) to be used for network topology characterization |46II50 |. 

Once the chemical system is assumed to be in steady-state, (C = 0): 


SV = 0 (142) 

where S is the stoichiometric matrix and V the specimens flux vector [mol.s — 1). The same problem can be 
derived for the mass-balance of each specimen: 


1 -1 
1 0 


-100 0 

-11-1 0 

-1 0 -1 -1 


^HCl 

’'U 2 

'"Cl 

"o 

"OH 

"H 2 O 


( 143 ) 


RV = 0 (144) 

where it can be shown that diag{kiCjCk) ■ S is equivalent to the 2nd term of eg 11411 being therefore an 
equivalent way of presenting reaction networks. If one considers the concentrations formulation, the reaction 
network dynamical system across the physical domain is given by: 
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N*Ndr2C + S I 


MdnCiC*U = 0 


or for a 1st order reaction kinetics: 


(145) 


N*NdflC + S ( 


■/ 


KdflCi = 0 


or by joining different mechamisms: 


(146) 


[ N^NdfiC + Si (g) [ MdnCiCfU + Sa® [ KdllCi = 0 (147) 

Jn Jn Jn 

Where Si and S 2 handle all the stoichiometric relationships between specimens. 

In chemical systems, network reconstruction is harder to cure when compared with biochemical data. In¬ 
formation is still scattered throughout publications and less efforts have been put into reconstructing chemical 
systems, such as in atmospheric science and foods. For example, Figure|4](b) presents part of known ascorbic acid 
(AA) degradation pathways |51[I52| . The full understanding of the AA degradation has major impact on both 
nutrition and quality of foods, but it still lacks the major mechanistic steps and thermodynamics. The same is 
valid for many important aging and degradation mechanism which involve oxidation IS3]. The reconstruction of 
this network implies the existence of high-throughput analytical chemistry dedicated facilities and bioinformatics, 
so that complex systems approaches can be applied to this research area [2]. 

As there is incomplete information, network simulation has to rely on flux analysis and measurements of flux 
rates instead of concentrations, kinetic rates, catalysis and Arrhenius activation energies. Considering that fluxes 
inside a triangular finite element is given by the shape function: 


v{Q,) = NiVi + NjVj + NkVk (148) 

where, Vi, Vj and Vk are the specimen flux at nodal positions i, j and fc; and the variational problem is resumed 
to: 


F(n) = ^ [ Sv-v{Q)dn = 
2 Jn 

-hi 


Ru dO, = 




{NiVi + NjVj + NkVkjdQ 


that once minimised for the node i, helds: 


Sv, 


* Jn 


{NiVi + NiNjVj + NiNkVk)dQ 


(149) 


(150) 


and performing for all nodal positions and chemical specimens, is possible to conclude the final matrix format: 


R® / N^NdfIV = 0 

Jn 

Where all stoichiometric relationships inside the FE space are respected, because: 


A = N^N = 


Nf NiNj NiNk 
NiNj Nf NjNk 
NiNk NjNk Nl 


(151) 


(152) 


and therefore, eg I151l in network of Figure 2] (a) is expanded to: 


r 

A 

A 

-A 

-A 

0 

0 

0 

/ 

0 

A 

0 

-A 

A 

A 

0 

Jn 

0 

A 

0 

-A 

0 

-A 

-A 


dfl ■ V = 0 


(153) 
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where 0 is a zero squared matrix, and V expands into a colum vector (21,1): 

[Vffi Vhj Vh^. ■ ■ ■ VH20i VH20fc] (154) 

Where all fluxes can be computed for any region of space. 

The same problem can be discretized using the stoichiometric matrix in eg 11411 where in complex chemical 
systems can be assembled from a knowledge base database table (Figure 14)1, where reactions, specimens, stoichio¬ 
metric factors, presence of catalysts, flux code and activation energies are cataloged, to obtain a linear system 
SV = 0, where in this example, V = [vri Vri Urs], and Vri = kiCHCnci, Vr 2 = k 2 CoCHCi and Vri = kiCoHCncu 
respectively. 

It can shown that inside any finite element, the set of equations became: 

S ® / N^NdflVe = 0 (155) 

Jq 

where for a triangular finite element, 


Vg — 

Being the solution for any given chemical network solve accross the physical domain as: 


L 


A 0 

A A 

0 -A 

-A -A 

0 A 

0 A 

0 0 


0 

A 

0 

-A 

0 

A 

A 


dO - V = 0 


(156) 


(157) 


In many cases, reaction mechanisms are not fully understood and incomplete analytical chemistry may not 
allow to derive all time-course dependencies in chemical systems. For engineering purposes, empirical pseudo¬ 
reaction steps can be assumed in incomplete reaction networks, such as for the ascorbic acid oxidation presented 
in Figure |4](b). The same formulation is possible to be presented to the pseudo-mechanistic network while there 
is not total knowledge about all reaction mechanisms (e.g. computational shelf-life dating |29|1. 


Effect of temperature and catalysts 

Pure chemical systems can be considered ’auto-regulated’ by thermodynamics, that is, mechanical properties, 
kinetic and equilibrium constants, activation energies and presence of catalysts. Chemical reactions dependence 
on temperature are generally modeled by the Arrhenius relationship: 


k = kj-ef X exp 



1 

T 



(158) 


where k and k^ef are the kinetic rates at temperature T and TrCf (K), respectively; Ea the Arrhenius 
activation energy {J.mol — 1, K — 1). The effect of catalysts can be reflected in the decrease of Ea, allowing the 
same reactions to occur at faster rates at lower temperatures. 

In order to reflect the effect of both temperature and catalysts, a weight matrix is possible to be deduced, as 
the fraction of the kinetic rate of a reaction step by it’s reference kinetic rate: 


Wi = 


ki 


demonstrating that under steady state the integration is given by S ■ diag(W) ■ V 
finite element formulation: 


(159) 


0 , with the corresponding 


S-diag(W)(gi [ N*NdfI-V = 0 (160) 

in 

where 0 < uii < -I- oo. Furthermore, when Wi = 0 the reaction step is deleted (e.g. deletion of a catalyst), 
Wi < 1 reactions are slower than the reference temperature, and Wi > 1 otherwise, enabling to study chemical 
systems under different environmental conditions. 
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'In-Silico’’ genome scale networks 

Modeling cellular growth had a significant impact on biotechnology in the pre-genomic era. Models with macro¬ 
scopic assumptions, also know as ’predictive microbiology’ (e.g. |54H56| 1 are still used due to their simplicity of 
assumptions and availability of information on kinetic data. FEM formulations were already derived for many of 
these models and can be found in [2]. 

The implementation of high-throughput methodologies in molecular biology (e.g. genome sequencing, elec¬ 
trophoresis, protein sequencing, mass spectroscopy, NMR), automated cellular manipulation (e.g. gene knock¬ 
out) [57] and the emergence of bioinformatics, provide that gene functions, protein specificity and partial metabolic 
networks are available in several species (e.g. ecoli, yeast and human) in databases such as, BioCyc [58], SGD |59| . 
KEGG [BD] , Reactome [BT], UniProt [B2]. With the increasing datasets, the development and update of holistic 
’in-silico’ genome-scale network draft models (GSM) has became possible to be automated |63II64 | for further 
validation by human experts to provide ’in-silico’ model organisms (Figure ??). 

There are three main types of ’in-silico’ GSM models: i) interaction network models; ii) steady-state stoichio¬ 
metric networks; and iii) dynamical models (e.g. ECELL [BB]). The latest are yet less used because of the lack of 
reliable ’in-vivo’ kinetic data, and therefore, interaction and steady state models are dominant in bioinformatics 
and systems biology analysis. Genome scale models can be further classified into non-compartmentalized and 
compartmentalized models (e.g. IND750 |30| . IMM940 [66]). The second class, accounts for metabolic networks 
contained in the different organelles and transport reactions between organelles, cytoplasm and extracellular 
space. Substantially complete models are available for ecoli (1260 genes, 2077 reactions, 690 of transport, 1039 
metabolites), s. cerevisiae (e.g.IND750, 750 genes, 648 metabolites, 1149 reactions, 297 of transport) and many 
other organisms in the BIGG database |67| . 


0.0.1 Flux-Balance Analysis 

Considering the example network inside an organism presented in Figure ??, at any given position of space inside 
a finite element domain, the concentration of metabolites of the ’in-silico’ organism can be given by: 


dt 


Sv + fiC = 0 


(161) 


where S is the stoichiometric matrix, v the metabolite flux (mol/s), the growth rate (s~^). In most 
conditions, as kinetic constants are not available ’in-vivo’, these models use the flux instead of the traditional 
kinetic constants. However, SkC = Sv, where v = kC at a given time or space position. Moreover, network 
studies assume pseudo steady-state conditions at a given time, that is, fluxes considered stable under short time 
periods, when compared to population growth and concentration of metabolites. The problem resumes to: 


which for the network model is: 


SV = 0 



VI 


• 0 • 


V 2 


0 


U3 


0 


U4 


0 


■^5 


0 




0 


^2 


0 


L ba J 


. 0 . 


(162) 


(163) 


Taking into consideration a consistent spacial gradient of the flux Va at any position, the solution is given by 
the minimization of the variational: 


V(H) = i (164) 

2 Jn 

and therefore, for a given metabolite i, the spacial solution is given by: 

S 0 / N^NdHVe = 0 (165) 

Jn 

where all network reactions are taken into account inside the hnite element space by using the Kronecker 
product with the stoichiometric matrix. Note that 14 is a column vector that spawns all vertices’s fluxes, such 
as, for a triangular finite element Ve = \vi,i • ■ • 64,i 64^ &4,)i]. 
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This simple formulation allows to perform FBA in conjunction with multi-physics FEM or CFD simulation 
in any biotechnological processes. In this sense, a state-of-the-art genome scale model analysis can be performed 
now with spatio-temporal resolution and in the complex scenario that modelers want to set-up for simulation, 
by integrating FEM solutions with systems biology to provide a genome scale diagnostic at any point of the EE 
mesh, such as the functionalities presented in )68| . 

State-of-the-art GSM were designed to assist molecular biology research, assuming chemostat conditions, 
and not for bioprocess or complex systems simulations. Today’s GSM cannot cope with: i) complex enzymatic 
mechanics; ii) kinetics and temperature effect; iii) dynamical states; iv) concentrations of metabolites; v) temporal 
and spacial resolution; vi) multi-physics phenomena are not taken into account (e.g. heat transfer, diffusion, 
fluid flow) and v) pathways are always assumed to be optimal, where control or thermodynamic restrictions are 
implemented by flux constrains. 

GSM provide today many applications in biotechnology, such as: i) flux balance analysis (FBA) for strain 
optimization; ii) network topological analysis, reliability, viability, structural homology; iii) derivation of pheno¬ 
type spaces for the exploration of biodiversity and biotechnological potential (Figure [5}. As these models do not 
hold a particular solution, both null space, convex analysis and optimization methods are applied to explore the 
solution space in chemostat conditions (e.g. MOMA, ROOM, genetic algorithms) [SS]. Furthermore, as solutions 
may converge into different regions of the phenotype, being necessary to develop new space basis, such as, the 
development of elementary flux analysis |69H71 | and extreme pathways |72H74| . 

The integration of GSM with FEM allows to overcome many of the previously mentioned barriers, allowing to 
perform genome-scale analysis of cells in the context of spatio-temporal conditions in a multi-physics environment 
[2]. Figure [5] exemplihes the integration. GSM are a set of incidence matrices, computationally derived from 
databases and cured with publications and expert analysis, relating genes to enzymes, enzymes and reactions, and, 
reactions to metabolites which given the stoichiometric relationships can be expressed as internal and boundary 
fluxes of metabolites. 

When deriving the GSM inside the finite element, the ’in-silico’ organism becomes dependent on the external 
conditions of nutrients, temperature, fluid flow, as well as, being affected by neighboring cells in any part of the 
physical and time domains. FEM considers that GSM is continuously discretized across the physical domain; at 
any point of the physical domain all metabolite fluxes and phenotype space is possible to be characterized, such 
as, for example the coordinates inside the convex hull given by the extreme pathways (see Figure [5] with limitless 
applications in biotechnology. 

Compartimented models 

In fully compartmentalized GSM models, each cellular organelle has an internal metabolic network, enzymes 
and associated genes. Gommon metabolites among compartments are linked by transport fluxes [SHIES]- In this 
reasoning, steady state equations resume to: 


SV-fTb = 0 (166) 

where T is the transport incidence matrix and b the boundary fluxes. After concatenation of all organelles 
metabolism and transport equations, cellular state inside a FE space is given by: 

S(g) / N^NdflV + T® [ N'^Ndnb = 0 (167) 

Jn Jn 

Being by this equation characterizes ’in-silico’ compartimentalized organisms at any region of the finite element 
space fl. 

Pheno-metabolomics 

Pheno-metabolomics plays a major role in post-genomic biotechnology. The exploration of the phenotype and 
metabolic capacities of organisms with the aid of both high-throughput methods in conjunction with genome 
scale models and complex systems simulation tools lies at the heart of pheno-metabolomics bioinformatics. Our 
research center has an important biodiversity yeast biobank, with especial emphasis on Saccharomyces cerevisiae 
isolates, and has been working in the characterization of S. cerevisiae over the last decade of yeast from different 
ecological contexts and geographical origins for their phenotype potential [75][7H]. 
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The pheno-metabolome of species is highly diversified, but most particular solutions of GSM have been 
restricted to the validation of simple, controlled experimental conditions | 30| which do not reflect the complexity 
of real-world bioprocess and natural conditions where yeasts evolved, lacking the design of new tools to both 
detect and derive new mechanisms as well as to cope with the dynamical complexity of cells. The use of GSM has 
been restricted to the assessment of the phenotype space derived from the stoichiometric matrix, being necessary 
to develop new approaches to fully explore the biodiversity of biobanks, evolution and adaptation mechanisms, as 
well as, the discovery of unknown mechanism by integration of GSM with both high-throughput signal processing, 
statistical computing, process analytical technology and computational simulation in order to be possible do derive 
the most correct definition of the phenotype space. 

One of the first approaches to define the phenotype of species was proposed as a non-negative linear com¬ 
bination of all relationships present in the stoichiometric matrix, holding all non-negative possible solutions of 
Sv = 0, when all fluxes Vi > 0 [77]. Such geometry is defined by the non-negative combination of a new vector 
basis, forming a convex hull defined by extreme rays (or pathways): 


71 

p = {w : w = ^ WiPi, Wi > 0} (168) 

i 

where p is the convex space (see Figure [T]) delimited by the extreme pathways Pi and Wi the coordinates 
projected into each Pi. Note that P is not an orthogonal basis, and only delimits the solution space of Sv = 0, 
being the vectors Pi presented in Figure[T]in the natural basis of Vi, which is not a practical visualization method 
once most GSM are hyper-dimensional. P can be obtained by the methodology presented in and hold 

important properties for the interpretation of the phenotype space: i) primary metabolism linked to boundary 
fluxes; ii) futile cycles with link to boundary fluxes; and iii) internal cycles. 

Inside a FE, the convex hull coordinates w of any point are possible to be described by the element shape 
function (or in any other basis): 


We = NiWi + NjWj + NkWk (169) 

allowing to apply finite element analysis (FEA) techniques do diagnose space differentiation in phenotype and 
metabolic state on the extreme pathways vector basis Pi. 


Spacio-temporal analysis 

Spacio-temporal analysis is perhaps one of the major advantages of joining FEM and GSM, becaming possible 
to analyze how the metabolic state evolves throughout space-time, as well as, to access how different phenotypes 
respond to different environment conditions. Previous sections already presented how to include fluxes (vi) 
and pheno-metabolome coordinates (wi) on a finite element domain. Such allows to analyze emergent patterns 
in cell communities and perform systems biology analysis | 68| at each region of space the cause of phenotype 
differences. Such tool will become more and more important, as cellular morphology may became manageable 
inside bioreactors iZHIlZSj. 

For instance, the use of the FEM allows to derive space vector gradients of both fluxes and phenotypes: 


dUe 

dt 


dNi 


-Ui -I- ■ 


dNi 


Z'Uj + 


dNk 


rUk 


(170) 


d{x,y,z) d{x,y,z) d{x,y,z) 
where Ue is the gradient property to be analysed across the physical domain. 

FEA may be used to further explore the phenotype dynamics, where for example the space derivate allows 
to determine geometrical changes in phenotypes across the FE domain (e.g. change rate and acceleration 


dUe dui duj duk 

—JT = + 

dt dt dt dt 


d^Ue 
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■ + Nk 


d Uk 
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(171) 

(172) 


Allowing to explore dynamically the molecular biology of different phenotypes, such as, the determination 
of the most important pathways and cellular functions at different stages, understand enzyme efficiency and 









27 


metabolic rates, regulation mechanisms and transcription rates in different contexts of cellular growth, as well as, 
understanding accelerations in phenotype changes or metabolic states as adaptations to changes in the environ¬ 
ment. Figure [5] resumes the use of the phenotype coordinates with FEM. 

As the solution of GSM equations is in many cases stochastic |68| . it is also important to be able to visualize 
the statistics of predictions in the FE domain. For example, is possible to derive both expected phenotype Ws 
and corresponding variance a'^{w) on a surface: 


f Nwds 

We = — - 
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a^(w) 


J^(Nw — We)^ds 
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(174) 


where s is the finite element surface area (m^) [2]. 

As this new approach may provide many possible solutions in the phenotype space, and therefore inverse 
FEM methods coupled with real-time and high-throughput experimental methods in molecular biology will be 
necessary to fine tune the numerical results of FEA analysis. Table [2] presents analogies between FEM-GSM 
and biological implications. Moreover, as dynamical results can be complex in terms of interpretation, pattern 
recognition recurring to compressed space coordinates may be more appropriate than direct visualization of fluxes 
and phenotype coordinates. 

The integration of FEM with reaction networks and genome scale networks will play an important role in 
the simulation and diagnostic of complex biological systems in the near future. Systems biology and systems 
chemistry lacked the possibility of integrating systems knowledge with multi-physics and multi-scale physics 
with 4D discretization that may enable in the future the computational assessment of phenotype tests, such 
as diagnostic the metabolic states under different growth media, emergence effects of gene deletion and stress 
factors, as well as bioengineering issues such as, reactor temperature, must composition and bioreactor design. 
This kind of tools will also open new possibilities in deriving and exploring the phenotype space for effective 
exploration of biobanks, providing critical informations for the decision of strain selection or improvement for 
a given biotechnological process. This manuscript is an introduction to the endless possibilities that are open 
for both study of complexity by FEM and network models and use of this methodology for the exploration of 
phenotypes, diagnosis, modeling, simulation and control of complex bioprocesses. 
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Figure 1. Key integration steps of reaction network models and finite elements: spacio-temporal 
discretization on FE space, time-course computation and results analysis in the phenotype space and 
fluxes, given different GSM configurations and environmental conditions. 
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Figure 2. Competitive inhibition inside a 3D finite element: reaction rates are fnnction of local 
concentration of specimens and temperatnre. 
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Figure 3. Example illustrating the FE concept with competitive inhibition: reaction rates are a 
continuous probabilistic function inside the FE space as function of concentration, given by 
approximation by the element shape function, providing a piecewise solution in the physical domain. 
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Figure 4. Chemical system reaction network: (a) mechanistic network representation of the set of 
chemical reactions in section 6; and (b) pseudo-mechanistic reaction network of ascorbic acid 
degradation in foods (adapted from 
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Figure 5. Main steps for the implementation of GSM in finite elements: i) development of draft 
models and human curation recurring to bibliography and experimental data; ii) development of the 
knowledge base for the ’in-silico’ organism implementation of the stoichiometry and transport matrices, 
control mechanisms and flux constrains; iii) assembling and solving FEM matrices for time-space 
resolution; iv) solving and analyzing results both in the phenotype space and physical domain imaging. 
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Figure 6. FEM applications for yeast pheno-metabolome exploration in biotechnology. 
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Table 1. Finite elements interpolation polynomials and shape functions 


Element 

Interpolation Polynomial 

Shape Function 

Linear Beam 

u{x) =71+ J2X 

?t = A^i??! + N2U2 

Linear Triangle 

u{x,y) =7i +722:+ 73?/ 

u = Niui + N2U2 + N3U3 

Linear Quadrilateral 

u(x, y) = 7i + 72X + 73?/ + 74a:?/ 

U = Nim + N2U2 + iV3?t3 + iV4??4 

Linear Tetrahedron 

u{x, y, z) = 71 + 72a; + 73?/ + 74Z 

U = NiUi + iV2??2 + N3U3 + iV4??4 

Linear Cube 

u{x, y, z) = 71 + 72X + 73?/ + 74Z 
+75X?/ + jexz + 'y^yz + ^sxyz 

u = Niui + ... + Nsug 

Quadratic Beam 

u{x) =71+ 72a; + 73X^ 

u = Nim + N2U2 

Quadratic Triangle 

u{x,y) =71 +722;+ 73?/ 

+742;^ + 73?/^ + 73a:?/ 

u = Niui + ... + Neug 

Quadratic Quadrilateral 

u{x, y) = 71 + 722; + 73?/ 

+742;^ + 75?/^ + 762;?/ 

+772;^y + 782;y^ 

u = Niui + ... + 

Quadratic Tetrahedron 

u{x, y, z) = 71 + 72a; + 73y + 742: 

+752;^ + 76y^ + 17 Z^ 

U = NiUi + . . . + NyU 7 

Quadratic cube 

u{x, y, z) = 71 + 72a; + 73y + 74Z 
+752;y + 752:2 + iryz + 782;y2 
+792;^ + 7ioy^ + 7112^ + 7i22;^y 
+713x^2 + 7i4xy^ + 7i5y^2 + li&xz^y 
+7i7y2^ + lisx'^yz + l\<s,xy^z + 7i62;y2^ 

u = Nim + ... + N20U20 


Table 2. Finite elements interpolation polynomials and shape functions 


Finite Element 

We 

We,cr(w) 

dWe 


dt 

d'^We 


dt"^ 

Vu> 


Biology _ 

Phenotype spacial distribution 
Phenotype statistical distribution 
Rate of cellular differentiation 
Rate of cellular adaptation 
Phenotype spacial differentiation vector 








